---
title: "outcome_ownstate_cov"
output: 
  html_document:
    toc: true
    toc_float:
      collapsed: true
      smooth_scroll: true
    depth: 3
    theme: cosmo
    highlight: pygments
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
```

```{r, eval = TRUE}
library(readxl)
ownstate <- read_excel("~/Documents/missouri_project/data/proportion_ownstate_guntrace.xlsx")
```

```{r, eval = TRUE}
library(dplyr)
library(car)
```

```{r, eval = TRUE}
# select MO and the eight bor states
data_ownstate <- ownstate %>% 
  filter(state=="Missouri" | state== "Iowa" | state=="Tennessee" | state=="Illinois" 
         | state=="Kentucky" | state=="Arkansas" | state=="Oklahoma" | state=="Kansas"
         | state=="Nebraska")
```

```{r, eval = TRUE}
# Remove the 2007 data
data_ownstate <- data_ownstate %>% filter(year!=2007)
```

```{r, eval = TRUE}
# Create an indicator variable about whether the subject is Missouri with year >= 2008
data_ownstate$MO_2008ge <- 0
data_ownstate$MO_2008ge[which(data_ownstate$year>=2008 & data_ownstate$state=="Missouri")] <- 1
```

```{r, eval = TRUE}
# Create another data frame with year only 2006, 2008-2013 data
data_ownstate2 <- data_ownstate %>% filter(year<=2013)
```

```{r, eval = TRUE}
# Turn state and year into factor variables
data_ownstate$year <- as.factor(data_ownstate$year)
data_ownstate$state <- as.factor(data_ownstate$state)

data_ownstate2$year <- as.factor(data_ownstate2$year)
data_ownstate2$state <- as.factor(data_ownstate2$state)
```

```{r, eval = TRUE}
### fit a linear model for year 2006, 2008-2019 data
fit1 <- lm(proportion~state+year+MO_2008ge-1,data_ownstate)
summary(fit1)
Anova(fit1)
```

```{r, eval = TRUE}
library(clubSandwich)
vcovCR(fit1, cluster = data_ownstate$state, type = "CR2")
```

```{r, eval = TRUE}
### fit a linear model for year 2006, 2008-2013 data
fit2 <- lm(proportion~state+year+MO_2008ge-1,data_ownstate2)
summary(fit2)
Anova(fit2)
```

```
vcovCR(fit2, cluster = data_ownstate2$state, type = "CR2")
```
